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(57) Abstract: A method for finding optimal filter co- 
efficients includes selecting a wavelet basis k parame- 
ters (102). The waveleet basis is reparameterized into 
k/2 rotation parameters (108) and factorized into a prod- 
uct of rotation and delay matrices (110). A data trans- 
form matrix is computed based on the product of ro- 
tation and delay matrices. The input data sequence is 
transformed into transformed data by the data transform 
matrix (1 12). The Jacobian of the data transform matrix 
and the input data sequence is determined and multiplied 
by the gradient vector with respect to the transformed 
data of an objective function. This product is compared 
to a predetermined criterium and if the predetermined 
criterium is not satisfied, a new set k/2 parameters val- 
ues are provided and the gradient descent is continued 
until the optimal k/2 parameters are found. The optimal 
filter coefficients are then calculated based on the opti- 
mal k/2 parameters (128). 
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TITLE OF THE INVENTION 
SIGNAL ADAPTIVE FILTER BANK OPTIMIZATION 

5 CROSS REFERENCE TO RELATED APPLICATIONS 

This application claims priority under 35 U.S.C. §19 (e) to 
US Provisional Patent Application Serial No. 60/269,678 filed 
February 20, 2001 the disclosure of which is incorporated by 
reference. 

10 

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR 

DEVELOPMENT 

15 

BACKGROUND OF THE INVENTION 
Filter coefficients can be optimized according to many 
different criterium. For example stop band attenuation, coding 
gain, degree of smoothness can be used as the criterium to 

20 optimize a filter. However, solving for these optimal filter 
coefficients requires solving a time domain problem that "is a 
non-linear constrained optimization problem that can be 
difficult to solve. In order to optimize these coefficients 
requires solving the non-linear constrained optimization 

25 equations iteratively to arrive at the optimum solutions. 
However, this iterative process is computationally expensive due 
to the difficulty of solving the equations. In addition, if the 
constraints are not satisfied then the solution is not 
invertible . 

30 Therefore, it would be advantageous to provide a method for 

solving for the optimal parameters for a filter that is not 
constrained and is more efficient computationally than the time 
domain problem. 

* 
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BRIEF SUMMARY OF THE INVENTION 

A method for finding optimal filter coefficients for a 
filter given an input data sequence and an objective function is 
disclosed. The method includes selecting a wavelet basis having 
k parameters and minimizing the predetermined objective function 
with respect to the k parameters. The wavelet basis is 
reparameterized into k/2 rotation parameters and factorized into 
a product of rotation and delay matrices. The k/2 rotation 
parameters are provided for the rotation matrices and a data 
transform matrix is computed based on the product of the 
rotation and delay matrices. The input data sequence is 
converted into transformed data by applying the data transform 
matrix to the input data. The Jacobian of the data transform 
matrix and the input data sequence is determined and multiplied 
by the gradient vector with respect to the transformed data of 
the objective function. This product is compared to a 
predetermined criterium and if the predetermined criterium is 
not satisfied, a new set of k/2 parameter values are provided 
and the gradient descent is continued until the optimal k/2 
parameters are found. The optimal filter coefficients are then 
calculated based on the optimal k/2 parameters. The wavelet 
basis may be selected from a wavelet packet library containing 
orthonormal wavelet packet bases, and in which the selected 
wavelet packet basis minimizes a given cost function, which can 
be an entropy function. 

A method for determining filter coefficients to form a 
filter to filter data of length N is disclosed, wherein the 
filter coefficients are optimized to minimize an objective 

■ 

function that measures a predetermined quality of the signal 
data. The method includes the steps of providing the number of 
coefficients, k, in the filter and selecting a wavelet packet 
basis. An objective function is provided and the first set of 
k/2 parameter values are also provided. A data transform matrix 
is formed as a function of the selected wavelet packet basis and 
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the k/2 parameter values. The data is transformed by 
multiplying the data transform matrix with the data, and the 
value of the objective function is calculated using the 

transformed data. The optimal set of values for the k/2 

.... ■ • 

5 parameters are then found based on the value of the objective 
function. 

In one embodiment, finding the optimal set of values for 
the k/2 parameters includes first determining if the objective 
function satisfies a first criteria. If the first criteria is 

10 satisfied, the optimal filter coefficients are calculated using 
the current set of k/2 parameters. If the objective function 
does not satisfy the first criteria a gradient steepest descent 
method is used to modify the k/2 parameters to a local minimum 
of the .objective function. The gradient steepest descent method 

15 is one in which the Jacobian of the data transform matrix is 
calculated and multiplied by the gradient of the value of the 
objective function with respect to the transformed data. If the 
product of the Jacobian and the gradient satisfies a 
predetermined criterium, the iterative process is stopped and 

20 the optimal filter coefficients are calculated based on the 
current set of k/2 parameters. If the product of the Jacobian 
and the gradient does not satisfy the predetermined criterium 
the values of the k/2 parameters are updated based on the value 
of the gradient, and the process is continued in an iterative 

25 manner until the criterium is satisfied and the optimal 
parameters found . 

Other forms, features and aspects of the above-described 
methods and system are described in the detailed description that 
follows. 

30 
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BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWING 
The invention will be more fully understood by reference to 
the following Detailed Description of the Invention in 

I ♦ ft • ■ w *» 

5 conjunction with the Drawing of which: 

Figs. 1A and IB together form a flowchart depicting one 
method of determining optimal filter parameters according to the 
presently disclosed invention; 

Fig. 2 is a diagram depicting a wavelet packet library of 
10 orthonormal wavelet packets; 

£ig. • 3 is a block diagram depicting a.' Given 1 s rotation 
block; and 

Fig. 4 is a block diagram depicting a cascade of lattice 
filters suitable for use with the method depicted in Fig. 1 

15 

DETAILED DESCRIPTION OF THE INVENTION 
A method for determining the optimal filter coefficients to 
form a filter of length M for filtering data of length N is 
disclosed. The filter is optimized to achieve an optimal value 

20 for an objective function that measures a predetermined, quality 
of the signal data. 

Fig. 1 depicts a flow chart for determining the optimal 
filter coefficients, wherein the number of filter coefficients, 
K, is selected, as depicted in step 102. The value of K can be 

25 selected based on the filter requirements, the system 
requirements, or a priori knowledge of the system response, or a 
combination of one or more these criterium. An objective 
function is provided, as depicted in step 103. The wavelet 
packet basis, i.e., the decomposition scheme is selected, as 

30 depicted in step 104, and a set of K coefficients are selected 
forming a predetermined mother wavelet, as depicted in step 106. 
The K coefficients are reparameterized into k/2. parameters that 
correspond to the mother wavelet in the lattice decomposition, 
as depicted in step 108. The data transform matrix is formed 
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using the decomposition scheme and the reparameterized k/2 
parameters, as depicted in step 109. The data transform matrix 
is factorized into a product of rotation and delay matrices, as 
depicted in step 110. The k/2 parameters are applied to the 
5 rotation matrices and used to provide the numerical values of 
the data transform matrix, as depicted in step 111. The data 
transform matrix is applied to a data sequence producing 
transformed data, as depicted in step 112. The objective 
function of the transformed data is calculated, as depicted in 

10 step 114, and compared to a first criteria, as depicted in step 
116...,: In^ the event that the calculated objective function 
satisfies the first criteria, control passes to step 128 and the 
optimal filter coefficients are calculated using the optimal 
parameters. In the event that the first criteria is not 

15 satisfied control passes to step 118. If the algorithm does not 
converge, a control loop can be used to limit the number of 
iterations and a new set of k/2 parameters can be provided and 
the process restarted. 

In step 118, The Jacobian of the f actorized-reparameterized 

20 second filter basis is calculated. The gradient of the 
objective function of the transformed data with respect to the 
transformed data is calculated, as depicted in step 120. The 
Jacobian and the gradient are then multiplied together to form 
the gradient of the objective function with respect to the set 

25 of k/2 parameters, as depicted in step 122. The gradient of the 
objective function with respect to the set of k/2 parameters is 
compared to a second criteria, as depicted in step 124. In the 
event that the gradient satisfies the second criteria, control 
passes to step 128 and the optimal filter coefficients are 

30 calculated from the first set of filter parameters. In the 
event that the gradient does not satisfy the criterium, control 

« 

passes to step 126. In step 126 the value of the k/2 parameters 

are updated as a function of the value of the gradient value and 

control passes to step 111, wherein the intervening steps are 
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repeated with the new set of filter parameters. This process is 
continued until either the objective function satisfies the 
first criteria or the gradient satisfies the second criteria. 

Alternatively, the initial values of the k/2 parameters can 

5 be determined other than by the time domain methods described 
above. For example, after executing steps 102, 103, and 104, 
the values of the k/2 parameters can be determined using other 
techniques. In one embodiment, the values of the k/2 parameters 
can be determined randomly and used in the optimization process. 

10 Alternatively, a user having a priori information or experience 
can select the values of the k/2 parameters accordingly. 
Alternatively, step 106 can also be executed and a mother 
wavelet selected and a set of a priori values for the k/2 
parameters can be supplied. In all of these alternative 

15 embodiments the optimization of the parameters can continue as 
described above with respect to steps 111-126. Step 128 is 
executed in the case where a digital filter architecture is 
used. In the case where the optimization and filtering is to be 
accomplished using lattice filters, there is no need to execute 

20 step 128 to find the coefficients of the digital filters. 

Filters can be optimized according to different criterium. 
The different criterium can be described according to a 
predetermined objective function. For example, stop-band 

attenuation, coding gain, and degree of smoothness are three 

25 possible criterium. Alternatively, the predetermined objective 
function can be based on morphological aspects of the signal to be 
filtered. For instance in medical monitoring applications, an 
objective function may be used to differentiate between normal 
breathing patterns and labored or wheezing breathing patterns, or 

30 between a normal cardiac rhythm and an abnormal cardiac rhythm. 
The objective function is application specific and developed 
according to a set of promulgated system requirements. 

In. general, a discrete wavelet transform includes low-pass 
and high-pass filters applied to a signal and the output is 
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down- sampled by a factor of two. The high-pass coefficients are 
retained and the high-pass and low-pass filters are applied to 
the low-pass coefficients until the length of the residual 
signal is smaller than the length of the filter. Accordingly , 

• » I t 

in general the maximum number of times a filter of length M+l, 
followed or preceded by a decimation by a factor of 2, can be 
applied to a data signal of length N+l, where N+l is a power of 
two, is given by: 



It should be noted that the data signal does not have to be a 
power of two for the optimal filter coefficients to be found. 
Although, the algorithm will operate more efficiently and faster 
for signal lengths that are a power of two it is not required. 
If the signal length is not a power of two, any of the well 
known techniques to extend the data set may be used. For 
example, zero padding where zeros are added to the end of the 
data signal to increase the signal length without interfering 
with the spectral analysis is an option as is the use of 
boundary wavelets to analyze the edges of the data signal. 
Other techniques that are known in signal processing may be used 
depending on the system requirements, the length of the data 
stream. to be analyzed and other design requirements. 

Consider a discrete wavelet transform operating on a signal 
with length that is a power of 2. As an illustration a filter 
of length 4 is operating on a signal of length 8. Let 

x=[x(0),x(1) 9 ..jc(7)] t be the original signal and c(0), c(l), c(2), and 
c(3) and d(0), d(l), d(2), d(3) be the low pass and high pass 
filter coefficients respectively. After applying the two 
filters to the original data, the result is down-sampled by two: 
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(3) 



the <y' s, the low pass 



(4) 



(5) 



These two steps can be combined and the transformed signal y is 
given by: 

y = Oc = C 2 C l x . (6) 

If Ci is orthonormal, then C 2 and C are also orthonormal. 
The conditions for orthonormality are given by: 

c(0) 2 +c(l) 2 +c(2) 2 +c(3) 2 =1 
c(0)c(2) + cQ)c(3) = 0 

d(0) 2 +d(l) 2 + d(2) 2 + d(3) 2 =1 
d(0) d(2) + d(l)d(3) = 0 

c(0)d(Q) + cQ)d(l) + c(2)d(2) + c(3)d(3) = 0 

c(0) d(2) + cQ)d(3) = 0 ( 9 ) 

c(2) d(0) + c(3)dQ) = 0 



(7) 



(8) 
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If the c's satisfy (7), then equations (8) and (9) can be solved 
using: 

d(k) = (-1) * c(3 -*),* = 0,...,3 . (10) 

As discussed above, the number of times filters of length M+l 

* 

5 can be applied to a signal of length N+l, where N+l is a power 
of two, is given by equation (1). Accordingly, the matrix C can 
be the product of Q orthonormal matrices: 

C = CqCq_ 1 ...C 2 C 1 , (11) 

where Q<Qmax and Qmax is given in equation (1) . 
10 The above example is intended for illustrative purposes 

only. In general, the constraints on the c's are given by: 

N 

£ c(n)c(n - 2k) = 5(k\ k = 0,. .., (N - 1) / 2 , (12) 

«=2Jfe 

the constraints on the d's are given by 
2 d(n)d(n - 2k) = S(k) 9 k = 0,..., (N - 1) / 2 , (13) 

n=2k 

15 and the relations between the c's and the d's are given by: 
Z!L c (" W w " 2 *> = <W> k = 0 - 1) / 2 

(14) 

EL* c (" - 2 *)<*(") = * = - 1) / 2 

* 

and where 5(k) is the Kronecker delta function. As above, when 
equation (12) is satisfied, equations (13) and (14) may be 
solved using: 

20 d(k) = (-1)* c(N -k),k = 0,..., N . ( 15 ) 

This leads to a constrained optimization problem that is 
difficult to solve for some applications. The coefficients 
c(0), c(l), ... c(N) are reparameterized so that the constraints 
in equation (12) are automatically satisfied. For example, 
25 equation (7) implies that: 

[c(0) + c(2)] 2 + [c(l) + c(3)] 2 = 1 , ( 16 ) 

which is automatically satisfied by setting 
c(0) + c(2) = cos(6\+0 2 ) 

c(l) + c(3) = sm(0 i +0 2 ) ' ( 
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Where the number of parameters have been reduced by one-half, 
from four to two in the illustrative example. Using the 
trigonometric formulas 

cos(ar + /?) = cos(a)cos(/?) - sin(a) sin(/?) 
sin(or + /?) = sin(a) cos(/?) + cos(a) sin(/?) 

leads to solving for the c's using: 



(18) 



c(O) = cos(0,)cos(0 2 ) 
c(l) = cos(0,)sin(0 2 ) 
c(2) = -sin(<9,)sin(<9 2 ) 
c(3) = sin(0,)cos(0 2 ) 



(19) 



For a length N filter, the orthonormality conditions imply that 



2>(«) 
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jnodd 



=1, 



(20) 



so that the c's may be reparameterized in terms of trigonometric 

* 

functions given by: 



(AM)/ 2 A*M)/2 \ 

£ C (2«) = cos j>„ 

/(AT-l)/2 \ 



z 

n=0 



c(2n + 1) = sin 



(21) 



2X 

V «=o j 



As can be seen, the right-hand side contains (N-l)/2 terms, 
but the expansions of the left-hand sides using generalized 
trigonometric addition formula contains 2 (N " 1)/2 terms. Although 
the reparameterized form of the filter uses less terms, 
distributing the trigonometric monomials to the various filter 
coefficients is difficult and complex. 

Distributing the trigonometric monomials to the various 
filter .coefficients is accomplished using lattice factorization 
of filter banks. A polyphase matrix of any two channel 
orthogonal filter bank can be factorized as: 



Ef = p(^)A(z)p(^ 2 ) • • • A(z)p(0 K ) 
where p(0)eO(2) and 



(22) 
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A(z) = 



1 0 

0 z' 1 

cos(0) sin(0) 

-sin(0) cos(0) 



(23) 



This procedure leads to a factorization of a wavelet 
transform, ' This factorization process is illustrated below for 
a filter having six (6) coefficients operating on a signal of 
length eight (8). Although illustrated for this example, the 
process can be extended for filters having a larger number of 
coefficients operating on signals having a longer length. 

Equation (22) can be rewritten as: 

H™=H?- x) A(z)p(0 K ). (24) 

As an illustrative example, a six coefficient filter using three 
angles will be provided. In particular, the polyphase matrices 
are given by: 



p p 



Id? d?l 



Hi = 



p 



(2) + z-'cf 



c< 3 > + z~ x cf + z~ 2 cf cf> + z-'cf + z~ 2 cf 
df> + z~ x df + z*d? d™ + z-'d? + z~ 2 d? 



(25) 



where the orthonormality is assured for H\, and the remaining 



polyphase matrices derived from H\ when 



cos(#, ) sin(#, ) 
- sin(#, ) cos(0, ) 



(26) 



Thus from equation (24), H 2 p is given by 

Hl=H l Mz) P {e 2 ) f 



(27) 



which can be written as: 
rff +z-V< 2) d^+z~ x df\ 



cos((9, ) sin(0, )" 
-sin(^j) cos(#,) 



1 0 

0 z- 1 



cos(# 2 ) sin(# 2 ) 
- sin(# 2 ) cos(# 2 )_ 



(28) 



Multiplying the right-hand side and solving for the value of the 
various coefficients on the left-hand side can be determined 
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using the reparameterized coefficients directly by equating 
powers of z. This leads to: 



cos(0! ) cos(<9 2 ) = c{° cos(0 2 ) 



c™ = cosfa )sinfe ) = c<'> sinfe ) 
cf = - sinfo )sin(0 2 ) = -<f sinfo ) 
cf > = sinfe )cos(0 2 )= c, (1) cos(0 2 ) 

and similarly for the d's: 
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■ sin(0, ) cos(0 2 ) = djp cos(& 2 ) 
•sin(<9 1 )siafe)=^ 1) sin(^ 2 ) 

■ cos(<9, )sin(0 2 ) = sin(0 2 ) 
«s(^j )cos(0 2 ) = d® cos(0 2 ) 



(30) 



cos(0, )cos(0 2 ) = rf™ cos(0 2 ) 
These relationships can be re-written as 
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(31) 



(32) 



for the relationships in equation (30) . As noted above in 
equation (24), the next polyphase matrix Hf* can be found as a 

function of . Multiplying the matrices and equating like 

powers of z as above yields: 
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(33) 
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Substituting equation (31) into equation (33) yields 
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= A 3 A 2 A^ 
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B T can be 



(34) 



B T =A*AlAl . Taking the transpose of equation (34) and combining 

equation (34) with a similar matrix equation for solving for the 
d's, yields: 



r P) r O) r 0) -(3) 

°5 °4 c 3 c 2 C l 

rf 5 (3) rf» </< 3) df> 

cos(0 3 ) -sin(0 3 ) 0 

0 0 sin(0 3 ) 

0 .0 cos(0 3 ) 

0 0 0 



,0) 



sin^ ) cos(^ ) 
cos(^ ) - sin(^ ) 

0 0 0 

cos(0 3 ) 0 0 
-sin(0 3 ) 0 0 

0 sin(0 3 ) cos(0 3 ) 



cos(0 2 ) 
0 



-sin(0 2 ) 0 0 
0 sin(0 2 ) cos(0 2 ) 



Thus, 



(35) 

function of the 



each coefficient can be found as a 
reparameterized angle coefficients. Performing the matrix 
multiplication in equation (35) defines each coefficient as a 
linear combination of two or more trigonometric functions. 
Reformulating equation (4) from above such that the c's and d' s 
alternate rows, the matrix Ci can be written as the product of 
three matrices: 
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C 5 




C 4 


c 3 


C 2 






c o 


0 


0 








0 


0 


0 


0 i 


0 


0 


d 5 


1 




d 3 


d 2 






do 


0 


0 








0 


0 


0 


0 i 


0 


0 


0 




0 








C 3 


c 2 








0 


0 


St 




0 


0 0 


0 


0 




0 




d 4 








*x 


d 0 





0 


0 






0 


0 1 


0 


0 






C 4 


0 


0 






C 2 




Co 




0 


0 


0 


0 


Si 


K x 0 


0 




1 


1 


0 


0 




^3 


m 

d 2 




d 0 




0 


0 


0 


0 




-S, 0 


0 


C 5 




C 4 








0 


0 


c i 






0 


0 


0 


0 


0 


0 s l 






1 


d 4 


* 

^3 






0 


0 




do. 




0 

IB 


0 


0 


0 


0 


0 K x 


-s 


"0 


1 


0 


0 


0 


0 


0 


0 






K 2 


0 


0 


0 


0 


0 


0 






0 


0 


1 


0 


0 


0 


0 


0 




*2 




*2 


0 


0 


0 


0 


0 


0 






0 


0 


0 


1 


0 


0 


0 


0 




0 


0 


^2 




0 


0 


0 


0 






0 


0 


p 


0 


1 


0 


0 


0 


* 


0 


0 


K 2 


-s 2 


0 


0 


0 


0 






0 


0 


0 


0 


0 


1 


0 


0 




0 




0 


0 


0 


*2 




0 


0 


* 




0 


0 


0 


0 


0 


0 


1 


0 




0 


0 


0 


0 


K 2 




0 


0 






0 


0 


0 


0 


0 


0 


0 


1 




0 


0 


0 


0 


0 


0 


- 


K 2 






1 


0 


0 


0 


0 


0 


0 


0_ 




_ 0 


0 


0 


0 


0 


0 


■^2 


Si. 






"0 


1 


0 


0 


0 


0 


0 


0" 






K 3 


0 


0 


0 


0 


0 


0 " 






0 


0 


1 


0 


0 


0 


0 


0 








*3 


0 


0 


0 


0 


0 


0 






0 


0 


0 


1 


0 


0 


0 


0 




0 


0 


s 3 


K 3 


0 


0 


0 


0 






0 


0 


0 


0 


1 


0 


0 


0 


* 


0 


0 




-s 3 


0 


0 


0 


0 






0 


0 


0 


0 


0 


1 


0 


0 




0 


0 


0 


0 




K 3 


0 


0 






0 


0 


0 


0 


0 


0 


1 


0 




0 


0 


0 


0 


*3 


-s 3 


0 


0 






0 


0 


0 


0 


0 


0 


0 


1 




0 


0 


0 


0 


0 


0 


*3 


K 3 






1 


0 


0 


0 


0 


0 


0 


0_ 




_ 0 


0 


0 


0 


0 


0 


K 3 


-s 3 _ 







(36) 

Where S k =sin(8 k ), K k =cos (9 k ) and all coefficients are assumed to 
be for the third polyphase matrix so that the superscripts have 
5 been dropped without any loss of generality. The matrix 
equation given in equation (36) can be expressed in the form of: 
C x =R(e x )SR(e 2 )SR(0 3 ) (37) 

where R(0) are wavelet transform matrices and S are shifting 
matrices. A matrix E is then applied to Ci to .separate the high- 
10 pass coefficients from the low-pass ones. 

The above is for illustrative purposes only, and signals 
with lengths of 16, 32, 64, or greater powers of two may be 
filtered. As an example, the left-hand side of equation (36) 
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10 



20 



25 



can be extended to filter a signal of length thirty-two (32). 
In this example, the left-hand side of equation (36), which is Ci 
in equation (11), becomes a 32x32 matrix in which the 
coefficients are shifted and wrapped around across all thirty- 
two columns. C 2 is formed in a block matrix form as follows: 



C 2 = 



^1(16x16) ^16x16 



L ^16*16 



16x16 _ 



(38) 



where d is the left-hand side 8x8 matrix of equation (36) 
extended to a 16x16 matrix by shifting and wrapping the filter 
coefficients as described above. In the illustrative example, 
Qmax is given by equation (1) and is equal to 3, C 3 is then 
given in block matrix form as: 



'1(8x8) 



0 



8x8 



0 



0 



16x16 



8x8 

0 



8x8 



16x16 



16x16 _ 



(39) 



where Ci is the left-hand side 8x8 matrix ,in equation (36) 
Thus, a wavelet transform of data of length 32 can be given by: 



15 C = C 3 C 2 C X = 



1(8x8) 



0 



8x8 



0 



0 



16x16 



8x8 

0 



8x8 



16x16 



16x16 



Q(16xl6) ^16x16 



_ ^16x16 



16x16 _ 



[A(32x32) ] ' 



(40) 



(41) 



This can also be expressed as: 
C^E Q R Q (0 l )S Q R e (0 2 y.S^ . 

As can be observed, the filter is a traditional straight 
wavelet transform. As is known, the straight wavelet transform 
continuously applies a low and a high pass filter to the data, 
down-samples the results by two (2), retaining the high pass 
results, then applies the high and low pass filters to the low 
pass results only. This is repeated until the filter size is 
too small to filter the remaining data. Alternatively, a subset 
of wavelet transforms that are less than the maximum number of 
transforms calculated in equation (1) may also be used to form 
the basis. 

Other wavelet bases may be used, however, and the possible 
combinations of the placement of the identity matrix and C x 
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matrix within each sub-block can be used. Fig. 2 provides a 
wavelet packet tree corresponding to a wavelet transform of 
input data 202 having a length of eight (8) . As discussed 
below, each subspace of the original signal space is filtered 

5 using the same predetermined high and low pass filters. The 
form of the predetermined high and low pass filters is selected 
depending upon the particular wavelet that is being used. The 
input data 202 is filtered by high and low pass filters, G and H 
respectively, into a first subspace having first and second 

10 orthonormal bases 204 and 206 respectively each having a length 
of four (4) . The first and second orthonormal bases 204 and 206 
are filtered by high and low pass filters, G and H respectively, 
forming a second subspace having four orthonormal bases 208, 
210, '212, and 214 each having a length of two (2). Each of the 

15 four orthonormal bases 208-214 are filtered by a high and low 
pass filter, G and H respectively, forming a third subspace 
having eight (8) orthonormal bases 216, 218, 220, 222, 224, 226, 
228, and 230 each having a length of one (1) . Each possible 
combination of non-overlapping subspaces of the original signal 

20 space forms a unique basis entry in the wavelet packet basis 
library, and in addition, each of the possible combinations is 
also an orthonormal basis. It is possible to identify the 
optimal orthonormal basis within the wavelet packet basis 
library that will provide for an optimal wavelet basis for 

25 representing that particular signal. 

A particularly efficient method of selecting the best basis 
is an entropy-based selection process. In this method, the 
entropy of a family of wavelet packet bases operating on a 
signal vector is calculated and the best basis is selected as 

30 the basis providing the minimum entropy. A suitable entropy 
function can be developed in which the sequences of data 
operated on by each basis in the wavelet packet library is 
compared by their respective rate of decay, i.e., the rate at 
which the elements within the sequence becomes negligible when 
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arranged in decreasing order. This allows the dimension of the 
signal to be computed as: 

d = exp 



~Y*P* l0 S2 Pn 

\ n J 



(42) 



where H x «f IMf and -^P n ^og 2 p n represents the entropy, or 

n 

5 information cost, of the processed signal. This leads to a 
proposition in which if two sequences {x n } and {x^} are compared 

so that a function corresponding function {p n } and \p n } are 

monotonically decreasing and if Y p» >V p[ then d ^ d' . 

Although entropy is one measure of concentration or efficiency 
10 of an expression, other cost functions are also possible 
depending upon the application and the need to discriminate and 
choose between special functions. 

The wavelet packet basis selected, which may or may not be 
the optimal basis for the particular signal, is applied to the 
sequence of data forming a transformed sequence of data. The 
particular parameters are optimized to select the optimal set of 
parameters and hence, the optimal set - of filter coefficients. 
In the preferred embodiment, a gradient descent can be used to 
optimize the parameters. In particular, the gradient of the 
20 objective function with respect to the parameter set can be 
calculated as: 

V^ = /JV^ (43) 

where J is the Jacobian matrix for Cx, § is the objective 
function, 9 is the set of parameters, and y is the input data 
25 sequence transformed by the selected basis. The Jacobian matrix 
is of the form: 

Je=[(deiC)x (^C)* (d*C)x] . (44) 

An explicit form of the derivatives with respect to the new 

parameters can be obtained using equation (41) : 

30 
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d j c ( 0 i 0 k ) = ipj c q Fq-\ -C l + C Q (dCg^ )..C t + ...C Q C q _ x ...(pjC x )= 

ER Q (0, )S Q ...d.R Q {ffj (0 K )+ER Q (0 l )S Q ..JjRq.JPj (0 K )-f 
• • • ER Q (0 l )S Q ...d J Rj \0j ) . JSxRx (0 K ) = 



er q(&i)Sq-R q 0j +^ \..S l R l (0 K )+ER e (0 1 )S Q ..J{Q_ l 



f 



2) 



.:£ v R l (0 K )+ 



J 2 



•••^i^i(^Ar) 



(45) 



where 



10 



15 



20 



dj=——, j=l,2,...,K and sin' (9)=cos (9) =sin (9+7i/2) , and 
00, 



cos' (9)=sin(9)=cos (9+7t/2) . It is possible to express the above 
derivatives in terms of rotations of the original angles, thus 
making the storage and computation of the derivatives more 
efficient. The derivative of a rotation block is: 



p{0)=p 



0+- 



2) L 



o 1 
-1 o 



(46) 



so that 

djR^D^j) 



(47) 



where Di is a block diagonal matrix. Thus, equation (45) 
becomes : 

djC(0 lt ...0 K ) = ER Q (0, )S Q ..JD Q R Q (0j Y^R, (0 K )+ 

ER Q (0, )S Q ..Dq^Rq., (0j (0 k )+ (48) 
ER e (0 l )S e ..JD l R l {ej\.S l R l (0 K ) 

Accordingly, with the appropriate choice of the matrices Ci, C2, 
. . . , C Q the above expression allows for the computation of the 
gradient with respect to the angles for any basis in the wavelet 
packet library. 

Thus, the Jacobian can be computed and stored efficiently 
and separately from the objective function gradient that is 
calculated. Accordingly, the Jacobian is efficiently computed 
for each set of parameters that are used in the gradient descent 
method. For each of the sets of variables that are used, the 
Jacobian can be calculated and stored according to equation 
(48) . 
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The gradient of the objective function (J> is calculated with 
respect to the transformed data sequence, i.e., the resultant 
data sequence from applying the selected wavelet basis to the 
input data sequence. This result is multiplied by the Jacobian 
5 matrix and the product of the Jacobian matrix and the gradient 
vector forms the gradient vector of the objective function with 
respect to the angle parameters. The gradient vector of the 
objective function with respect to the angle parameters provides 
a vector that is the steepest descent from the starting point to 
10 a local minimum. If the gradient vector is equal to zero, or 
less than a gradient threshold level, a local minimum has been 
reached, or nearly reached, and the objective function is 
calculated for the current angle parameters. If the objective 
function is less than a predetermined objective function 
15 criterium, the current angle parameters are the optimal angle 
parameters and the filter coefficients can be calculated as 
described above. If the gradient vector of the objective 
function with respect to the angle parameters does not indicate 
that it is at or near a local minimum, or if the objective 
20 function exceeds the predetermined objective function threshold, 
then a predetermined adjustment is made to the current set of 
angle parameters. This process is repeated until an optimal set 
of angle parameters is found. The actual value of the gradient 
threshold level, the predetermined objective function threshold, 
and the predetermined adjustment to the angle parameters are 
dependent on the particular system requirements. 

* 

Fig. 3 depicts a block diagram of a rotation block suitable 
for rotating a pair of input signals according to the rotation: 

(49) 

30 where Ke=cos (0) and Se=sin(0) . The rotation block depicted in Fig. 
3 is considered to have a rotation angle of G such that the input 
signals will be rotated through an angle of 0. The rotation block 
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300 includes first and second inputs oti and a 2 . Input ai is 
multiplied by the cos (0i) and sin (61) , and input a 2 is multiplied 
by the cos (9i) and -sinOx). The products aiCos(Oi) and -a 2 sin(0i) 
are added together in summation module 302 to provide output ai R . 
5 Similarly, the products <x 2 cos (0i) and otisinOi) .are added together 
in summation module 304 to provide output <x 2 R . The rotation block 
depicted in Fig. 3, in combination with delay modules and down- 
sample modules, forms the basic building block of the lattice 

filters . When used in a wavelet system, these lattice filters can 

* 

10 form the high and low pass filters that are used to determine the 

■ 

wavelet transform of the input data. 

Fig. 4 depicts one form of a cascade lattice filter that is 
suitable for use with a filter of length four (4) using two 
decomposition steps. The architecture and layout of the various 

15 sub-systems is based on the straight wavelet transform. In this 

1 

formulation the length of the signal can be arbitrary since the 
values of the signal are continuously modified from left to right. 

Fig. 4 is a circuit equivalent of steps 111-118 of the flow 
chart depicted in Fig. 1. The lattice filter depicted in Fig. 4 

20 is of length four (4) for illustrative purposes only. The filter 
bank depicted therein can be modified to include more parameters, 
i.e., angles, more decomposition steps, i.e., larger Q, and 
different bases from the wavelet packet table by adjusting the 
overall architecture in accordance with the system requirements. 

25 In addition, as the angles are updated during the optimization 
process described with respect to Fig. 1, the corresponding 
rotation angles in the various lattice filters will be updated 
accordingly. 

In particular, Fig. 4 depicts a wavelet transform/ Jacobian 
30 module 400 that receives a sequence of data on line 401 and 
filters the received data with filter module 402. Filter module 
402 is a two stage lattice filter that includes two rotation 
modules as depicted in Fig. 3 that have rotation angles equal to 
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* 

81 and 0 2 respectively. Filter module 402 provides a first high 
pass data output and a first low pass data output, 403 and 405 
respectively, that represent the high pass and low pass filtered 
input data respectively. The low pass data output 405 is 
subsequently filtered by filter module 404 that includes two 
rotation j blocks that have rotation angles of 9i and 9 2 
respectively. The filter module 404 provides a second high pass 
data output and a second low pass data output 407 and 409 
respectively. The first high pass data output 403 and the second 
high pass data output 407 and the second low pass data output 409 
together form the transformed data output values y(0), y(l), and 
y ( 2 ) respectively . 

The Jacobian can also be computed simultaneously with the 
transformed data vector y(0), y(l), and y(2). In particular, the 
data on line 413 is negated by module 414 and provided as a first 
input to provided module 412. The data on line 411 is provided as 
the second input to module 412. Module 412 in combination with 
the portion 402A of filter module 402 form a lattice filter module 
having first and second rotation blocks having rotation angles of 
0i and 0 2 respectively. Module 412 provides a first data output 
415 and a second data output 417, wherein the second data output 
417 is subsequently filtered by lattice filter module 410 that 
provides a third and fourth data output 419 and 421. 

The high pass filter data output 403 is provided to lattice 
filter module 408 that has two rotation modules that have rotation 

angles of 0i and 0 2 respectively. Filter module 408 filters the 
high pass data output 403 and provides a fifth data output 423 and 
a sixth data output 425. The data on line 427 is provided to 
module 406, which in combination with the first portion 404A of 
lattice filter module 404 also forms a lattice filter having two 

rotation blocks that have rotation angles of Q x and 0 2 
respectively. Module 406 provides a seventh data output 429 and 
an eighth data output 431. 

* 
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The Jacobian values are formed as linear combinations of the 
various output data. In particular, J(0) is formed as the sum of 
the low pass data filter output 405 that is negated by module 416 
and the transformed output Y(0) combined in summation module 420. 
5 J(l) is formed as the sum of the fifth data output line 423 and 
the transformed output value Y (2) , negated by module 418, combined 
in the summation module 422. J(2) is formed as the sum of the 
sixth data output line 425 and the transformed output value Y(l) 
combined in the summation module 424. J (3) is formed as the sum 
of the first data output line 415 and the transformed output value 
Y(0) combined in the summation module 426. J (4) is formed as the 
sum of the third data output line 419 and the seventh data output 
on line 429 combined in the summation module 428. J(5) is formed 
as the sum of the fourth data output line 421 and the eighth 
15 output data line 431 combined in the summation module 430. the 
data that is provided on the outputs J(0) , J(l) J(5) corresponds 
to the type of decomposition matrix that is used. For example, in 
the straight wavelet transform depicted in equations (38) -(40) the 
various Y output values correspond to the various blocks of the 
20 block diagonal matrix. In this example using 64 data input 
signals, J(0) provides 32 values, J(l) provides 16 values,. J(2) 
provides 8 values, J (3) provides 4 values, J (4) provides 2 values, 
and J(l) provides 1 value. As can be seen, the values provided by 
each output corresponds to the position of the identity matrix in 
25 the straight wavelet transform block matrices. 

It is convenient to require that the mean of the high pass 
filter used herein be zero to ensure that no DC component passes 
through the filter without significant attenuation. In the time 
domain formulation a zero mean high pass filter satisfies the zero 
30 mean filter criteria: 

S(-l) B c(i\T-n) = 0. (50) 



n=0 



where c is the coefficient of the high pass filter as discussed 
above. This condition when translated into the lattice angle 
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formulation and combined with the inequalities in equation (21) 
becomes : 



£(-!)" C (W-«)= t l c(2n)-f i c(2n + l) = 



COS 



H=0 



H=0 



-sin 



=o. 



(51) 



The above is true if 



(52) 



where k is equal to N/2, and N is the number of time domain filter 

Thus, the search for the optimal filter 



coefficients. 



coefficients becomes a constrained optimization problem that is 
expressed as: 



9 . = 

rain 



argmin 



(57) 



The gradient of objective function can easily be expressed on this 
plane, where the objective function is given by: 



(58) 



V 4 M J 

then the gradient of the objective function provided in equation 
(58) is given by: 



1 



0 -1 



0 ... 1 -1 



\ 



\ 



* M J 



(59) 



As can be seen, constraining the high pass filter to be a zero 
mean filter has the advantage of reducing the complexity of the 
optimization problem by one from k to k-1 degrees of freedom. 

Those of ordinary skill in the art should further appreciate 
that variations to and modification of the above-described methods 
for providing optimal filter coefficients may be made without 
departing from the inventive concepts disclosed herein. 
Accordingly, the invention should be viewed as limited solely by 
the scope spirit of the appended claims. 
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What is claimed is: 

1. A method for determining filter coefficients to form a 
filter to filter data of length N, wherein the filter 
coefficients are optimized to minimize an objective function 
that measures a predetermined quality of the signal data, the 
method comprising the steps of: 

providing a number of coefficients, K, in the filter; 

selecting a wavelet packet basis; 

providing an objective function; 
s providing a first set of k/2 parameter values; 

forming a data transform matrix as a function of wavelet 
packet basis and the k/2 parameter values; 

calculating transformed data by multiplying the data 
transform matrix with the signal data; 

calculating the value of the objective function based on 
the transformed data; and 

finding the optimal set of values for the k/2 parameters. 

2. The method of claim 1 wherein the step of finding the 
optimal set of values for the k/2 parameters includes: 

determining if the value of the objective function 
satisfies a first criteria; 

in the event that the value of the objective function 
satisfies the first criteria, setting the optimal set of the 
values of the k/2 parameters equal to the current value of the 

» 

first set of k/2 parameters; 

in the event that the value of the objective function does 
not satisfy the first criteria, 

calculating the Jacobian of the data transform matrix 

based on the data and the current values of the k/2 

parameters; 

calculating the gradient of the calculated objective 
function with respect to the transformed input data; 
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multiplying the Jacobian with the gradient of the 
calculated objective function with respect to the transformed 
input data to form the gradient of the parameters; 

determining if the gradient of the parameters 
5 satisfies a second threshold value; 

in the event that the gradient of the parameters satisfies 
the second criteria, calculating the optimal filter coefficients 
based on the current values of the k/2 parameters; 

in the event that gradient does not satisfy the second 
10 threshold value, forming a new set of k/2 parameter values by 
updating the value of the current value of the- k/2 parameters 

« 

delta value; 

returning to the forming a data transform matrix step 
and re-execute the intervening steps. 



15 



25 



3. The method of claim 1 wherein the step of selecting a 
wavelet packet basis includes selecting an optimal wavelet 

* 

basis . 



20 4. The method of claim 3 wherein the optimal wavelet basis is 
an orthonormal basis. 



5. The method of claim 4 wherein the optimal, .orthonormal basis 
is selected from a wavelet packet library. 



6. The method of claim 5 wherein the optimal orthonormal basis 
selected from a wavelet packet library is selected to minimize a 

* 

predetermined cost function. 



30 7. The method of claim 6 wherein the cost function is an 
entropy function of the wavelet packet basis applied to the 
input data sequence . 
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8. The method of claim 7 wherein the entropy function is equal 
to log 2 p n . 



n 



9. The method of claim 1 wherein the step of providing a first 
5 set of k/2 parameters includes providing k/2 lattice filter 

* 

angles . 



10 



10. The method of claim 9 wherein the reparameterized basis is 
a function of at least one trigonometric function. 



11. The method of claim 1 wherein the data transform matrix is 
of the form C = E Q R Q fa )S Q R Q (ff 2 )..S Q R Q fa \.JE X R X fa )S X R X fa \.S X R X fa ) with 



log. 



r M+\ \ 



15 



20 



12. The method of claim 2 wherein the form of the Jacobian 
matrix of the data transform matrix is of the form 

ER Q (0 X )S Q ...d.R Q (Oj \.S X R X (0 K ) + ER Q (0 X )S Q ...djR Q _ x (Oj )..S X R X (0 K ) + 
...JLRq (9 x )S Q ...djRj(0j)>..S x R x (0 K ) 

13. The method of claim 12 wherein the form of the Jacobian 
matrix of the data transform matrix is of the form 



* 



( A 
0J+- 



14. The method of claim 13 wherein the Jacobian matrix of the 

data transform matrix is of the form 

ER Q (9, )S Q .. £> Q R Q [dj (0 i )+ER Q (0 1 )S e ...D e _ l R Q _ 1 (fl,)..^ (<?,)+ 



25 
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15- The method of claim 1 wherein step (d) providing a first set 
of k/2 parameter values includes providing random values for the 
first set of k/2 parameters. 

5 16. The method of claim 1 wherein the step of providing a first 
set of k/2 parameter values includes using a priori information. 

17. The method of claim 1 wherein the step of providing a first 
set of k/2 parameter values includes the steps of: 

10 selecting a mother wavelet; 

selecting a set of K coefficients corresponding to the 
mother wavelet; and 

computing the k/2 parameters that correspond to the mother 
wavelet in a lattice decomposition. 

15 

18. The method of claim 1 wherein the step of providing a first 
set of k/2 parameter values includes the steps of: 

selecting a mother wavelet; 

selecting a set of K coefficients corresponding to the 
20 mother wavelet; and 

providing k/2 predetermined parameters that correspond to 
the mother wavelet in a lattice decomposition. 

19. The method of claim 1 wherein the step of providing a first 
25 set of k/2 parameter values includes providing the first set of 

k/2 parameter values that satisfy a zero mean filter criteria. 

r 

20. The method of claim 19 wherein the zero mean filter 

30 criteria is Y&, = — . 

M 4 

21. The method of claim 20 wherein the step of finding the 
optimal set of values for the k/2 parameters includes: 
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defining the objective function as 

( k-\ ^ 

providing the gradient of the objective 
respect to the k/2 parameters as: 



function with 



1 



0 



0 -1 



1 -1 



7t 



4 



k-i \ 

M J 



setting the gradient to a threshold value; and 
solving the gradient for a set of k/2 parameters that 
satisfy the gradient equation. 



-28- 



WO 02/071256 



PCT/US02/04420 



1/5 



Select the number of 
coefficients, K, In the filter 
102 
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Select an objective function 

103 


• 


r 



Select the wavelet packet 
basis, I.e., the decomposition 
scheme 
104 



Provide a first set of K 
coefficient values 
106 
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Reparameterize the K 
coefficients into K/2 
parameters 
108 
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Form the data transform 
matrix 
109 


I — i 





Factorize the data transform 

matrix into a product of 
rotations and delay matrices 

110 



Apply the K/2 parameters to 
the Factorlzed data transform 
matrix 
111 



0 



Apply the data transform 
matrix to the input data 
sequence forming 
transformed data 
112 



Calculate the objective 
function of the transformed 
data 
114 



Does the objective 
function satisfy the first 
criteria 
116 



Yes 

i 



Calculate the K optimal filter 
coefficients from the set of 
K/2 parameters 

128 •■ 
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Fig. 1A 
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0 

i 

Calculate the Jacobian of the 
data transform matrix with 
respect to the input data and 
the current set of K/2 - 
parameters 
118 



1 




Calculate the < 
objective functi 
to the trans! 
1i 


gradient of the 
on with respect 
farmed data 
JO 



Calculate the gradient of the objective 
function with respect to the current set 
of K/2 parameters as the product of the 
Jacobian and the gradient with respect 
to the transformed data 
122 



foes the gradient with respect t< 
the current set of K/2 parameters 
satisfy a second criteria 
124 



•Yes- 



Update the value of the K/2 
parameters based on the value 
of the gradient with respect to 
the current set of K/2 
parameters 
126 
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